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Abstract 



This paper describes the f ortran program rhad which performs a numerical eval- 
uation of the photon-induced hadronic _R-ratio, R{s), related to the cross section for 
electron-positron annihilation, for a given center-of-mass energy ^/s. In rhad the 
state-of-the-art perturbative corrections to R{s) are implemented and the running 
and decoupling of the strong coupling constant and the quark masses is automatically 
treated consistently. Several options allow for a flexible use of the program. 



PROGRAM SUMMARY 

Title of program: rhad 

Available from: 

http : // www . rhad . de/ 

Computer for which the program is designed and others on which it is operable: Any 
work-station or PC where f ortran is running. 

Operating system or monitor under which the program has been tested: Alpha, Linux, 
Solaris 

No. of bytes in distributed program including test data etc.: 168000 
Distribution format: ASCII 

Keywords: Hadronic i?-ratio, perturbative QCD, electron-positron annihilation, run- 
ning and decoupling of ag 

Nature of physical problem: The hadronic i?-ratio R{s) is a fundamental quantitiy in 
high energy physics. It is defined as the ratio of the inclusive cross section a{e~^e~ — ^ 
hadrons) and the point cross section apt = 47ra;^/(3s). It is well-defined both from 
the experimental and the theoretical side. R{s) belongs to the few physical quantities 
for which high-order perturbative calculations have been performed (partial results 
up to order exist!). Mass effects from real and virtual quarks, the evolution of 
the MS parameters, in particular in the presence of thresholds, and other subtleties 
lead to fairly complex results in high orders. Thus it is important to provide a 
comprehensive collection of formulas in order to make them available to non-experts. 

Method of solution: rhad is a compilation of all currently available perturbative 
QCD corrections to the quantity R{s). Several options are provided which allow for 
a flexible use. In addition, rhad contains routines which perform the running and 
decoupling of the strong coupling constant. Thus only the center-of-mass energy has 
to be provided in order to determine R{s). 

Restrictions on the complexity of the problem: The applicability of rhad is restricted 
to the perturbative energy regions and does not cover the narrow and broad reso- 
nances. 

Typical running time: The typical runtime is of the order of fractions of a second. 
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LONG WRITE-UP 



1 General structure of R{s) 

We are considering the fully inclusive production of quark pairs in e~^e~ annihilation 
(for a review see Ref. [1]). The tree-level diagram for this process in shown in Fig. 1 (a). 
At leading order (LO) perturbative QCD (pQCD), the cross section as a function of the 
squared center-of-mass (c.m.s.) energy, s, has thresholds at the points s = ^rrt^ with 
Q = d,u, s, c, b, t. However, close to threshold, fixed-order pQCD is no longer applicable. 
Below threshold, non-perturbative effects lead to the formation of bound states of the 
quark-anti-quark pair, which appear as more or less sharp peaks in the cross section.^ In 
our perturbative framework of charm and bottom production, these narrow resonances 
cannot be described, and arc therefore not included in rhad. In addition, we have to spare 
out the region between the physical threshold Sq" and the beginning of the more or less 
flat continuum region at Sq"', where the cross section exhibits rapid variations. In the 
case of charm production, for example, the limits would be i/s|f" ^ 2m d ~ 3.73 GeV 
and i/s^ 4.8 GeV, while for bottom production they are ^ 2mB « 10.52 GeV 

and Y^s^ ~ 11.2 GeV. Even though rhad provides a numerical value in these threshold 
regions, one should not try to give this value any physical significance; rhad will print a 
warning that reminds the user of this fact. 

Furthermore, we have to exclude the low-energy region, below about y^Smin = 1-8 GeV, 
where the validity of perturbation theory is doubtfiil. Thus, mass effects from u-, d- and 
s-quarks arc negligible, which is why we will consider these quarks as massless throughout 
the paper. 

Instead of the cross section cr(e+e~ — hadrons), we will express the results in terms of 
the so-called hadronic i?-ratio, 

, cr(e+e~ —> hadrons) 

R{s) = \ , 1 

where 

aie-e- = '-^ (2) 

is the high-energy limit for the photon-mediated lowest order muon pair production. In this 
paper we concentrate on the contributions induced by photon exchange. Z-boson exchange 

is only relevant at high energies. There, however, apart from the QCD corrections^, also the 
electro- weak corrections become important [2]; these will not be addressed in this paper. 
Nevertheless, rhad can be used to evaluate the vector contribution to the production of 
top quarks in the same way as for charm and bottom production. 

We write R{s) as a sum of contributions from individual quark flavors, plus a so-called 

^Due to the large width no bound state is formed in the top quark system. However, a peak in the 

cross section remains around s — 4m1. 

^For a discussion of corrections to the axial- vector correlator see Ref. [1]. 
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singlet term: 



Q=d,u,s,c,b,t 

The "non-singlet" contributions Rq{s) are all proportional to the square of the respective 
quark charge. The (numerically small) singlet piece -Rsing appears for the first time at 
order and will be discussed in Sect. 6. 

Rq{s) is defined to be zero below the continuum region of QQ production (s^q^), i.e., we 
do not attempt to describe the resonance regime: 

Rq{s<s'^'^) = 0. (4) 



Our main concern arc the radiative corrections to Rq{s), predominantly arising from 
gluons exchanged between, or emitted from, the QQ pair, but the lowest order QED 
effects, though very small, will also be included. Thus we write 



Rq{s) = e{s - s] 



thr\ 
Q ) 



n>0 



(5) 



This formula generally describes the QCD/QED effects to Rq{s), if SRq^^{s) contains 
all contributions that vanish when electro-magnetic corrections are ignored. Both ag and 
Rq\s) explicitely depend on the renormalization scale /i, while Rq{s) is invariant under 
/x-variation, up to and including the order of perturbation theory that is being considered. 

Since the QCD corrections only affect the quarks in Fig. 1 (a), the leptonic part of the 
diagrams can be factored out and one remains with corrections to the 'jQQ vertex, and 
real radiation of quarks and gluons. Beyond order as, the most convenient way to evaluate 
the fully inclusive rate is to compute the photon self energy and relate it to the rate 
7* QQ + X through the optical theorem, i.e., by taking the imaginary part in the 
time-like region of the external momentum. Explicitely, 



R{s) = 127rlmn(s + ie) , 



(6) 



where 

n(^^) = ^(-5,. + ^)n,.(,) (7) 

is the transversal part of the photon polarization function, n^iy(g). 

The different contributions to R{s), as we classify them in our paper, will be illustrated by 
sample diagrams contributing to n(g^) in what follows. The LO terms are thus represented 
by Fig. 1 (h) . Since the imaginary part of these diagrams is obtained by the sum of all 
possible cuts, it includes real and virtual contributions simultaneously. 

A final remark concerning the counting of loops is in order. We understand by the ri-loop 
QCD contribution to R{s) the sum of real and virtual corrections to order a". Computing 
R{s) via the optical theorem at n-loop level thus requires the evaluation of the imaginary 
part of the photon self-energy to n -|- 1 loops. 
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2 Running and decoupling of as and the quark masses 



As already mentioned in the previous section, the quantity R{s) depends on the strong 
couphng constant ag, the heavy quark masses mg, and the renormahzation scale /x. All 
of them arc specified in the subroutine parameters which is discussed in Appendix C. 
The dependence of R{s) on fj, is explicit in terms of logarithms, but also implicit as, for 
example, in the argument of ccg. At higher orders in perturbation theory, one needs to 
specify the scheme in which ag and the quark masses niQ are evaluated. For a^, we will 
always adopt the MS scheme which is commonly used in QCD calculations. Concerning 
the quark masses, it is generally more appropriate to use the pole mass in the energy 
region close to the production threshold of the quarks. On the other hand, the use of 
the running mass and setting 1x^ = 3 resums part of the potentially large logarithms 
Inrnq/s in the high-energy region. For this reason, rhad provides the logical parameter 
Imsbar that allows the user to switch between the pole and the MS mass definition for 
the evaluation of R{s). If not stated otherwise, we denote by mq a generic quark mass 
in what follows. In those cases where it is necessary to distinguish between the pole and 
MS quark mass, we denote them by Mq and ffiq = 'inQ^\iJ.), respectively. The conversion 
formula between Mq and mg is given in App. E.4, Eq. (67). The "scale-invariant" mass 
is defined recursively as fiiQifiiQ) = fhQ-^\fhQ^^), where the number of active flavors is 
equal to rij = 4 for charm, Uf = 5 for bottom, and Uf = 6 for top quarks. 

If one chooses to evaluate R{s) in terms of MS masses, the input required by rhad is 
a^^\Mz) and the scale invariant masses fhQ{fhQ) (Q = c,b,t). In addition, one needs 
to specify the renormahzation scale /x and the c.m.s. energy squared, s. From that, rhad 
will determine the number of active flavors nf, as well as the parameters a^^^\ii) and 
ffiQ^^^ji), which will be used for the evaluation of Ri^s). As already mentioned before, a 
natural choice (and the default value in rhad) for the renormahzation scale is /U^ = s. 

The number of active flavors is determined through the threshold variables s*'^^, s^^^ , and 
gthr -yy^hen s < s*^'', it assumes the smallest possible value, = 3, and it increases by 
one every time one of the threshold variables is crossed, up to = 6 for s > s*^"". 

On the other hand, in the definition of a^"^\fj,) one has the freedom to choose "matching 
scales" /Xc, fib, and /it, at which the transition from ri/ to ra/± 1 is performed. For example, 
assume that n/ = 4, and the renormahzation scale is set to some specific value n = p.- The 

procedure to compute ai^\fl) at n-loop order as implemented in rhad is as follows: In a 
first step the renormahzation group equation (RGE) for the strong coupling (here and in 
what follows, we refer to App. E.4 for the definition of the coefficients /3, Q, 7^, and C^), 

"'d?^ = ,3'»" («i"''w) = - g /?<"'' (^^l^j . (8) 

is solved numerically in order to obtain al (iJ'b) from the input value al {Mz)- The use 
of the decoupling relation 

af^-'\nb) = (C.)'«^\/x,) (9) 
to (n — l)-loop order, leads to a^f^{ii},). Finally, applying Eq. (8) a second time gives the 
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value for ai'^\fi). The described procedure is performed automatically in rhad using the 
subroutine rundecalpha (see App.D). 

The evaluation of rriQ^ (ji) from the input fhQ^rfiQ) proceeds completely analogously. The 
RGE that governs the running of the quark masses in the MS scheme reads: 



(10) 



Combining Eqs. (8) and (10) leads to 



N ^ c(q;,. (/i)/7r) 

rriQ {n) = mQ (/^o) ^^^^^^ , (11) 



with 



c{x) = x''" 1 1 + (ci - 6ico) 3^ + ^ [(ci - feico)^ + C2 - 6ici + 6f Co - 62C0] ic'' 
+ 



^(ci - bicof + ^(ci - 61 Co) (c2 - 61C1 + 6? Co - 62C0) 



+ ^ (c3 - ^'iC2 + &?ci - b2Ci - bfco + 2 5162C0 - 63C0) + 0{x'^) I 



(12) 



Cj = 



If R(s) is evaluated at n-loop order, the curly bracket of Eq. (12) has to be evaluated up 
to, and including, the term oc a;"~^. 

The matching between rfiQ^ ^ and ^ is determined through the equation 

m^^^-'^ = Cmfn^^^\ (13) 

where the function (j^ can again be found in App. E.4. The procedure for the evaluation 
of fnQ^\ii) from the input mQ^ffiQ) in rhad is called rundecmass (see App.D). 

Note that the evaluation of fnQ^\ii) is only required if the user decides to use MS masses 
for the evaluation of R{s). In this case, in parameters . f (see App. C), 

Imsbar = . true . 

and the values for masse, massb, masst have to be set to the scale invariant masses fhc{ific), 
rhbirhf,) , ■fhti'fht), respectively. If one chooses to use pole masses instead, one defines 

Imsbar = .false. 

and sets masse, massb, masst equal to the on-shell masses M^, M^, Mt, respectively, rhad 
will then simply use these values throughout the calculation. 
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There is one exception where we do not allow a choice between on-shcU and MS mass: 
For "power-suppressed" corrections, i.e., double-bubble diagrams that contain a heavy 
secondary quark loop (see Fig. 3 (e) with m2 S> -y/s, mi), the heavy quark mass is al- 
ways inserted in the on-shell scheme. If the input for the heavy mass is provided in the 
MS scheme, the corresponding on-shell value is evaluated within rhad by the subroutine 
iimis2nios (see App. D), using the proper conversion relations [3, 4, 5]. 

3 Notation and tree-level result 

(a) (b) 




Figure 1: (a) Tree-level graph for e+e — > QQ. (6) One-loop photon polarization function. 



(a) (6) 




Figure 2: Diagrams contributing to R{s) at one-loop level. 

The tree level result for Rq{s) is part of every introductory course in quantum field theory. 
It reads 

i?g)(s) = neV§rf (s,mQ), with r'^\s , Mq) = ^ (3 - v^) , (14) 

where ric = 3 is the number of colors, and Vq is the electric charge of the quark Q in units 
of the proton charge. I.e., 

Vu = Vc = Vt = +'^, Vd = V. = Vfc = -^. (15) 
VQ denotes the quark velocity, 

where ^/s is the c.m.s. energy. Note that we define vq in terms of the on-shell mass, Mq 
(see Sect. 2). This becomes relevant at higher orders, where we refer to the definition of 
Eq. (16). 
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4 One-loop result 



Sample diagrams that contribute to the one-loop approximation of R{s) are shown in 
Fig. 2. Their imaginary part has been evaluated a long time ago in the context of QED [6]. 
In its most compact form, the result reads 



R^^\s) = ncVl^r^^\s,mQ), 



(17) 



with 



r'i\s,MQ) = CF 



3~v 



+ SvQ In 



Q 



(1 + v^q) Li2 (/) + 2Li2 (p) - Inpln 



l-vl 33 + 22vl - 74 ^^l + VQ , l^VQ - 



VQ In VQ + 



8(3 - v^q) 



+ 



1-VQ 4(3 - V, 



(18) 



wherep = {^ — vq)I(\-\-vq) and Cp = 4/3. In the limit vq — 1, one obtains Vy = 3Cf/4 
which, after taking into account the coupling and color factors, leads to the following 
expression for the QED corrections: 



6R 



QED 



3 a 

^4 TT 



(19) 



The tree-level and one-loop functions and their correspondence in rhad are summarized 
in Tab. 1. 



notation 


in rhad 


diagrams 


(0) 


rvO 


Fig. 1 {h) 




rvl 


Fig. 2 



Table 1: Contributions to R{s) at tree-level and one-loop order. 



5 Two-loop result 

Starting from two-loop order, a general analytic expression for R(s) is no longer available. 
Typical diagrams are shown in Fig. 3. The full set can be obtained from these diagrams 
by attaching the two external photon lines to the quark lines at arbitrary points. Only 
certain classes of diagrams have been evaluated in closed form. Nevertheless, the analytic 
evaluation of different kinematical limits, combined with appropriate interpolation among 
these limits, have resulted in extremely accurate semi-analytical approximations to the 
full result. Thus, for all practical purposes, the full mass dependence is available at order 

Two main new features occur at the two- loop level that are absent at lower orders. First, 
there are contributions that are due to the non-Abelian character of QCD and have no 
correspondence in QED. These are diagrams with a three-gluon coupling (see Fig. 3 (c) and 
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(d); the four-gluon couphng occurs for the first time at three-loop order). The second new 
feature is that a second quark hne can appear, with the effect that the result may depend 
on two different quark masses mi and m2 (see Fig. 3 (e)). 

It is clear that the simplest set of diagrams to evaluate are the ones with self-energy 
insertions in a gluon propagator, Fig. 3 (d), (e). In fact, the cases with massless insertions 
(i.e., gluons and massless quarks) are known in closed analytical form, for general values 
of s and mi [7, 8]. The same is true for mi = 0, m2 7^ [9] (labels in accordance 
with Fig. 3(e)). The case mi = m2 is known in terms of a two-dimensional integral 
representation [7, 10]. For mi <C m2, the coefficient of the leading term in an expansion 
around s/m^ is known in closed form [11]. For m,i = [12], this term was shown 
to approximate the exact s/m^ dependence [9] extremely well, even up to the threshold 
s = 4m2. 

The diagrams without self-energy insertions on gluon lines have not yet been evaluated 
in closed form for general values of mq and s. The massless limit has been known for 

quite some time [13]. Subsequently, mass corrections in the high energy limit have been 
computed: The tuq/s terms can be obtained by a simple Taylor expansion of the in- 
tegrand [14]; the mq/s^ terms were derived from the MS operator product expansion 
of n(g^), combined with renormalization group relations [15]; and the evaluation of the 
higher order terms (up to rriq /s^) was achieved by systematically expanding the Feynman 
integrals [16]. 

However, one should note that the convergence of such a high-energy expansion is not 
guaranteed below the highest threshold of the diagrams under consideration. For example, 
the non-planar diagram in Fig. 3 (6) has a four-quark cut at = (4mQ)^, meaning that 
the expansion around ruq/s — formally converges only above s = {4mq)'^. Nevertheless, 
in practice one often observes also satisfactory convergence below the threshold [16]. 

A result that is valid in the full kinematic range was obtained in Ref. [17] (see also Ref. [18]). 
To this aim, the high-energy limit was combined with up to eight moments of the polariza- 
tion function obtained through an expansion for 0. This information was combined 
with additional input from the threshold region around q^ = 4mg, using a conformal 
mapping and Fade approximation in order to arrive at a three- loop expression for n(g^). 
For all known practical applications, this approximate result is equivalent to an analytic 
expression for Il{q'^) and, after taking the imaginary part, for R{s). 

Let us now explicitely parameterize the two-loop contribution to R{s). We can write 



it;g)(s) = neV^ 



rv{s, mq) + ry (s, mq) + V ry'(s, mq,mq] 



(20) 



where the sum runs over all quark flavors, ry oc C| denotes contributions from Abelian 
diagrams, ry oc CaCf comes from non-planar and non- Abelian diagrams, and ry' oc CpT 
("double-bubble") arises from diagrams with two separate quark lines. Cp = 4/3, Ca = 3 
and T = 1/2 are color factors of QCD. For later reference, it is convenient to introduce 
additional functions that refer to specific limits of ry": 

rv'is,mq,0) = ry{s,mq) , ry'{s,mq,mq) = ry{s,mq) . (21) 

See also Tab. 2 for more details. Note that ry" introduces contributions from secondary 
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quarks q other than both virtual and real. They arise through the splitting of gluons 
into a quark-anti-quark pair. 

In the massless case, one has the simple expressions 



where L^^ = hi{s/iJi^) and Cs ~ 1.20206. The expressions in the massive case are given in 



In the following, the expressions adopted for the individual contributions are described, 
and their origin is given. 

• For (s, mq) and (s, mq) we use the full mass dependence as given by the Pade 
results of Eqs. (36), (39) [17]. 

• For the double-bubble contribution ry'[s,mi,m2), we have to distinguish several 



— mi = m2 = mq: For y/s < ^mq, we use the analytic formula as given 
in Eq. (44) [7]. For y/s > ^mq, the full mass dependence is known in terms of 
a two-fold integral representation, see Eq. (46) [7]. For the sake of speed, how- 
ever, we use the high energy expansion of Eq. (48) [16]. Numerical differences 
to the integral representation are completely negligible. Nevertheless, the user 
may switch to the integral formula by setting Irvctexp = .false, in rhad. 

— m2 < mi: We apply an expansion for ^ mf (< s/4): 

1. At 0{m2), the full mi-dependence is used in the form of the analytic result 
of Eq. (41) (for x = L) [7]. 

2. At 0(771,2), we neglect the mi dependence. It turns out that this contribu- 
tion vanishes which can be seen in Eq. (52). 

3. Higher orders in m2 are neglected. 

- (2mi)^ < s < (2?TT,2)^: If mi = 0, we use the analytical result of Eq. (49) [9], 
including the full dependence on m2- If "'i ^ 0. wc use the leading term in 
1 / rr?! , keeping its full dependence, see Eq. (54) . This expression has been 
obtained in Rcf. [11]. 

- (2mi)2 < (2m2)^ < s: If mi = 0, we use the analytical result of Eq. (49) [9], 
including the full dependence on m2- If mi 7^ 0, we apply an expansion in the 
limit ml < m|(< s/4): 

1. At 0(m5), we use the analytical result of Eq. (49) [9] for the full m2 de- 
pendence. 

2. At 0{ml), we neglect all effects from non-zero m2- The corresponding 
expression for ry'(s,mi,0) is given in Eq. (53). 

3. Higher orders in mi are neglected. 




(22) 



App.E.l. 



cases: 
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Note that we set the ratios m^/m^ and m^/m^ to zero and consider mass corrections due 
to the presence of an additional Ught quark only for c-quark effects in 6-quark production. 

For the sake of completeness, let us remark that also the two-loop corrections of order 
aag and order for massless quarks [19] are known. Numerically, these contributions are 
very small. Nevertheless, we include the mixed corrections of order aag into rhad, which 
modifies Eq. (19) to 

6Rq _neVQ--^l--^J . (23) 

Higher order QED or mixed QED/QCD effects are numerically irrelevant and will be ne- 
glected [20]. 

(a) (6) (c) 




(d) (e) 




Figure 3: Classes of diagrams contributing to R{s) at two-loop level. (For the issue of 
counting loops, see the discussion at the end of Sect. 1.) 



notation 


in rhad 


diagrams 


mi 


m2 


4 


rvcf 


(«): 


'"Q 






rvca 


(6), (c), (d) 








rvdb 


(e) 


mi 


m2 


ry 


rvnl 


(e) 


mQ 





T 


rvct 


(e) 


rriQ 





Table 2: Different contributions to R{s) at two-loop order. The column "diagrams" refers 
to Fig. 3, where a sample diagram is shown for each class. 



6 Three-loop result 

At order a^, a new class of diagrams contributes to -Rhad, the so-called "singlet diagrams." 
They differ from the non-singlet diagrams in the sense that, in the singlet case, the external 
currents are not connected by a common quark line. A typical diagram is shown in 
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(a) 



id) 





Figure 4: Classes of diagrams contributing to R{s) at three-loop level. 



notation 


in rhad 


diagrams 


mi 


m2 


ms 




rv3 


(a), (6), (c) 


mg 


mQ or 


mQ or 




delr03 









mg or 


y,sinK 


rvSsing 


id) 


mi 


m2 





Table 3: Different contributions to R{s) at three-loop order (as far as available). The 
column "diagrams" refers to Fig. 4, where a sample diagram is shown for each contribution. 



Fig. 4 (ci). The non-singlet diagrams (see Fig. 4 (a)-(c)) are numerically dominant, but 
rhad includes both contributions. We will describe them in more detail in what follows. 



6.1 Non-singlet contributions 

The knowledge of R{s) at three-loop level is restricted to the high energy limit. The 
massless limit was obtained for the first time in Ref. [21] and later confirmed through an 
independent calculation in Ref. [22] . 

Mass corrections are known up to 0{mQ/s^): The quadratic terms were obtained from 

renormalization group identities [23], while the quartic mass terms were evaluated by 
applying the methods of Ref. [15] at three-loop level [24]. 

Let us briefly discuss how we treat mass effects from massive inner quark loops. At three- 
loop level, there can be three different quark types in one diagram, with, say mi > m2 > 
ms. The approximation that we apply is to take all terms including m\ and m\, neglecting 
mixed terms of order mfm^. Mass effects from ms are also neglected. In addition, we do 
not consider power-suppressed terms, i.e. diagrams that contain a quark q with s < s^^"^. 

With these approximations in mind, we can write the three-loop contribution as 



i?g)(s) = neV§ 



'•^y\s,mQ,nf) + 8rf{s,mq,nf) 

q+Q 



(24) 



{s,mq,nf) denotes the mass corrections due to massive inner quark loops different 
from Q. According to the approximations above, it vanishes for s < s^^^ and also for 
mq = 0. 

Note that n/ itself depends on the c.m.s. energy. Assume, for example, that s < s^^^. 
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Then we have Uf = 3 and 

i?§) (s) = ne rg) (s, 0, 3) , Q = n, d, s . 
For s*^"^ < s < s^^"^, we have n/ = 4 and thus the non-zero contributions are 



and (s) = nc V^^ r{^^ (s, m^, 4) . 

For sj;'"' < s < Sj*^'', it is = 5 and 



4\s,0A)+Sr^'"{s,mcA) 



(3), 



Q = u,d, s, 



(25) 



(26) 



i?§)(s)=n,V§ 



r|^^ {s, 0, 5) + 5r^^^ [s, rric, 5) + (5r^^^ (s, m^, 5) 



Q = u,d, s, 



(s) = V2 (s, me, 5) + 5ri'^ (s, m^, 5) 
iif ) is) = Tie Vl [s, mb, 5) + (5rf ^ (s, me, 5) 
Finally, for s > sf^'^, all masses except for mt are neglected, and, with n/ = 6, we have 



(27) 



i?g^(s) = neV^ 



s,0, 6) + (5rQ^^(s,mt,6) 



Q = u,d, c, s, b, 



Rl''{s) = ncVi4\s,mt,6). 



(28) 



As mentioned before, in Eq. (27) it is understood that all c-quark mass effects are included 
only at 0{ml). 

The expressions for ry\s,mQ,nf) and 6r^\s,mQ,nf) are listed in App. E.2. 



6.2 Singlet contributions 

As already mentioned above, singlet diagrams are characterized by the fact that each of 
the external photons couples to a separate quark line (see, e.g., Fig. 4((i)). The singlet 
contribution arises first at three-loop level; at two-loop level, where the two quark lines 
are connected by only two gluons, it is zero due to Furry's theorem. 



We write the singlet contribution to R{s) in the following way: 

-Rsing(s) =nc^Vq Vqi r{f^ing(s, rUq, Ulq') , 

q,q' 



(29) 



where the sum over q and q' runs over all active quark flavors. E.g., for -y/i = 12 GcV 
it includes all quarks except for the top quark. Again, only the high energy expansion 
is known. Quadratic mass terms are absent, and the expression including quartic mass 
terms reads: 

(3) , X 55 5 ^ m^ / 10 25 \ m^ / 10 25 ^ \ 

Note that, if m^ 7^ m2, we neglect the lighter of the two masses. Since this is the lowest 
order where singlet terms occur, Eq. (30) is the same in the MS and in the on-shell mass 
scheme. 
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7 Four- loop result 



The knowledge of the four-loop contribution to R{s) is very limited. Only contributions of 
diagrams with two and three quark-loop insertions have been calculated. Sample diagrams 
are shown in Fig. 5 (c) and (d). The mass effects are evaluated up to order mq/s, so that 
the results are strictly only valid in the high-energy limit. 



(a) 



id) 





Figure 5: Classes of diagrams contributing to R{s) at four-loop level. 



We write the four-loop contribution as 



with 



V 



P.(4),0L^ 



(s) + rif fS'^'^^s) + nj r^'^'^'^is) + n} r^'"''''{s) 



ncVQrv\s,mQ) , 



2 ^(4),2L, 



3 ^(4),3L, 



(31) 



(32) 



where the bar indicates that the expressions are estimates rather than approximations or 
even exact results, n/ is the number of active flavors, varying with the c.m.s. energy as 
discussed in Sect. 2. 

^(4),3L ^j^^ renormalon-type contribution with three massless quark insertions on the 
gluon propagator, see Fig. 5(fi). It has been evaluated in Rcf. [25], where the general 
structure of the terms of order as(asnj)" was derived. Recently, the analytic expression 

for rl^^'"^^ became available [26]. It required the evaluation of massless four-loop two-point 
functions in combination with the method of Ref. [22] (see also Ref. [18]) to derive the 
imaginary part of the five-loop contributions as shown in Fig. 5 (c) . The same method has 
been used in combination with the technique derived in Ref. [23] for the computation of 
the quadratic mass corrections ry2^^ and Vy}^^^ [27]. 

Estimates for the full massless four-loop result have been known before and are still a 
subject of interest [28, 29, 30]. We determine fl^\s) and rl)'^'^{s) by subtracting the 
known nj and nj contributions from the estimates of Ref. [28] atn/ = l,...,6, and fitting 
the resulting six "data points" by a linear function in nf. One finds, for /x^ = s, 



.{4),0L 



,{4),2L 



(s) = -1.86 • 10^ 



.(4),1L 



(s) = 21.3, 



-7.97 • 10" 



„(4),3L 



(s) = 2.15 • 10" 



(33) 



The logarithmic contributions follow from the lower order terms through renormalization 



group invariance and are collected in App. E.3. Once the exact results for Tq 
r^^'^^ become available, the approximate ones can easily be replaced in rhad. 



{4),0L 



and 
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Using the same method for the quadratic mass terms in combination with the estimates 
of Ref. [27] we obtain, in the MS scheme 

= 7.11 • 103 , ftl'^^s) = -1.43 ■ 10^ , 
4^)'^^.)=49.1, 4^)'=^^.) = -0.204. ^ ^ 

In App. E.3, the corresponding result for on-shell quark masses is hsted together with the 
logarithmic contributions. 

For completeness, the function implemented in rhad is listed in Tab. 4. Note that at 
0{ag), the corrections analoguous to Sr^^ are neglected. 



notation 


in rhad 




rv4 



Table 4: Contribution to R{s) at four-loop order. 



8 Evaluating R{s) 

In this section we describe how rhad evaluates R{s) in the individual energy regions. Let 
us first recall that defines the lowest value of the c.m.s. energy squared, s, at which 
the perturbative treatment of QQ production is allowed; Rq{s) is defined to be zero for 
s < Sq^^. The physical threshold for QQ production is at Sq^; the user is advised, however, 
to disregard the results of rhad in the region between Sq^ and Sq"' {Q = c, b, t) . 

The evaluation of the number of active flavors n/, the strong coupling constant ai"^ ■*(//), 

and, if required, the MS quark masses friQ^\ii), from the input values s, /x, ai'\Mz)^ and 
fhQ{fhQ), has been described in Sect. 2. Let us now look in more detail at the specific 
contributions that enter R{s) at certain values of s. 

For s < sj°™, i.e., below the production threshold for two Dq mesons, we have Uf = 3. 
The charm, bottom and top quark masses are decoupled and thus their contribution goes 
like a^s/niq {q = c, b, t), see Eqs. (49,54). All other correction terms are evaluated in the 
massless limit. 

In the region s^^^ < s < s|,°^ one has n/ = 4. ai^^ (/x) is evaluated from the input ai^^ {Mz), 

and, in the MS mass scheme, fh'^\fi) is obtained from the input quantity fhcifhc) (see 
Sect. 2). In the on-shell mass scheme, the input value Mc is used unchanged throughout 
the calculation. The full charm quark mass dependence (inasmuch it is known) is taken into 
account. Bottom and top quark are decoupled and enter through the power-suppressed 
terms, Eqs. (49,54), at order a^. 

At c.m.s. energies above and below Sj°^, the situation is similar to the previously 

described region, but with Uf = 5. Note, however, that in the MS scheme not only 
fh^^\fi) is needed for the numerical evaluation, but also fh'^\iJ,). It is obtained from the 
input fhc{fhc) using the routine rundecmass (see Sect. D). 
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= 2 GeV y/s = b GeV 











OIC16I 






U[S) 





0.1180 


2.0012 







0.1180 


1.2199 


3.2220 


1 


0.2726 


2.1747 




1 


0.2025 


1.4956 


3.6267 


2 


0.2981 


2.2223 




2 


0.2120 


1.5950 


3.7467 


3 


0.2984 


2.2049 




3 


0.2123 


1.5825 


3.7272 


4 


0.2973 


2.1835 




4 


0.2122 


1.5467 


3.6867 



Table 5: Hadronic i?-ratio at ^/s = 2 GeV (nj = 3, left) and ^/s = 5 GeV (nj = 4, right). 

Note the change in ag which is consistently evaluated from ai^\Mz)- The parameter 
settings are given in App. C.2. 



^/s = 12 GeV 



order 


as 


Rc{s) 


Rb{s) 


R{s) 





0.1180 


1.3304 


0.2675 


3.6001 


1 


0.1667 


1.4199 


0.3496 


3.8778 


2 


0.1706 


1.4269 


0.3803 


3.9265 


3 


0.1709 


1.4198 


0.3797 


3.9148 


4 


0.1709 


1.4158 


0.3756 


3.9050 



Table 6: Hadronic i?-ratio at ^/s = 12 GeV (nj = 5). 



For s > s\^^, a^^\n) and inf' {ii) arc needed for the evaluation of R{s). At these energies, 
no mass corrections from lighter quarks are considered. 

Tab. 5 and 6 show the results for R{s) with the default settings (see Sect. C), for different 
values of the c.m.s. energy. The individual rows correspond the different orders of per- 
turbation theory. In the case of ^/s = 5 GeV and ^/s = 12 GeV, we also show the partial 
contributions to R{s) arising from massive quarks. The second column contains the value 
of a^J^^\-\fs), with the appropriate n/. 

Tab. 7 contains the results for as{s), R{s), Rc{s) and Rb{s) for the energy region between 
y/s = 1.8 GeV and ^/s = 20 GeV, sparing out those parts where perturbation theory is 
not applicable. We adopted the default values of rhad given in App. C with iord = 4. In 
Fig. 6 the data for R{s) (full curve) are shown in graphical form. The shaded bands mark 
the uncertainty of R{s) induced by the variation of as{Mz), /x, Mc and Mf, in the range 



,(6) 



as{Mz) = 0.118 ±0.003, 

Mc = (1.65 ±0.15) GeV, 
Mb = (4.75 ± 0.2) GeV , 



(35) 
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(GeV) 




Rc{s) 


Rb{s) 


R{s) 


1.80 


0.3148 


0.0000 


0.0000 


2.1895 


2.20 


0.2833 


0.0000 


0.0000 


2.1780 


2.60 


0.2620 


0.0000 


0.0000 


2.1684 


3.00 


0.2463 


0.0000 


0.0000 


2.1G07 


3.40 


0.2342 


0.0000 


0.0000 


2.1544 


3.73 


0.2260 


0.0000 


0.0000 


2.1499 


4.80 


0.2150 


1.5681 


0.0000 


3.7097 


5.00 


0.2122 


1.5467 


0.0000 


3.6867 


6.00 


0.2007 


1.1812 


0.0000 


3.6174 


7.00 


0.1919 


1.4556 


0.0000 


3.5836 


8.00 


0.1850 


1.4399 


0.0000 


3.5637 


9.00 


0.1793 


1.4301 


0.0000 


3.5505 


10.00 


0.1744 


1.4235 


0.0000 


3.5410 


10.52 


0.1722 


1.4209 


0.0000 


3.5371 


11.20 


0.1736 


1.4186 


0.3793 


3.9133 


12.00 


0.1709 


1.4158 


0.3756 


3.9050 


13.00 


0.1679 


1.4129 


0.3717 


3.8964 


14.00 


0.1652 


1.4106 


0.3685 


3.8891 


15.00 


0.1628 


1.4087 


0.3659 


3.8830 


16.00 


0.1606 


1.4070 


0.3636 


3.8778 


17.00 


0.1586 


1.4056 


0.3618 


3.8732 


18.00 


0.1567 


1.4044 


0.3602 


3.8692 


19.00 


0.1550 


1.4033 


0.3589 


3.8657 


20.00 


0.1534 


1.4023 


0.3577 


3.8626 



Table 7: as{s), Rc{s), Rb{s), and R{s) for various values of s. 
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2 3 456789 10 20 

Vs"(GeV) 

Figure 6: R{s) for 1.8 GeV < -y/s < 20 GeV. The central curve (solid line) is obtained using 
the default values specified in App. C and corresponds to the results given in Tab. 7. The 
error bands are obtaind by varying as{Mz), /x, Mc and Mb within the limits of Eq. (35). 
Not shown are the perturbatively inaccessible regions between s^™ and s*^^ (Q = c, b, see 
Sect. 1). 



9 A typical program 

It is instructive to look at a typical program that evaluates R{s). We set the c.m.s. energy 
to ^/s = 12 GeV, and split the result according to the contributions from the individual 
quarks. The occuring functions will be described in more detail in App.B, but their 
functionality should be clear from the following program. 

Input. The corresponding f ortran program would look as follows. 

program example 

implicit real*8(a-h,m-z) 
implicit integer (i,j) 
implicit character*60(k) 
implicit logical (1) 

include ' common . f ' 
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sqrts = 12. dO 

scms = sqrts*sqrts 



call parameters (scms) 
call init(scms) 



rail 


= rhad(scins) 


ru = 


ruqrk(scms) 


rd = 


rdqrk(scms) 


rs = 


rsqrk(scms) 


rc = 


rcqrk(scms) 


rb = 


rbqrk(scms) 


rt = 


rtqrk(scms) 


rsg = 


= rsinglet (scms) 


rem = 


= rqed(scms) 



print* , 


'R_had = 


' ,rall 


print* , 


'R_u 


' ,ru 


print* , 


'R_d 


',rd 


print* , 


'R_s 


' ,rs 


print* , 


'R_c 


' ,rc 


print* , 


'R_b 


',rb 


print* , 


'R_t 


' ,rt 


print* , 


'R_sing = 


' .rsg 


print* , 


'R_QED = 


' ,rem 



end 



Output. If Iverbose is set to .true, the output reads: 



rhad.f — Version 1.00 
by Robert Harlander and Matthias Steinhauser 
December 2002 

Order of calculation: 4 

Scales : 

sqrt(s) = 12.000 GeV 

mu = 12.000 GeV 

thrc = 4.800 GeV 

thrb = 11.200 GeV 

thrt = 360.000 GeV 

Number of active flavors: 
nf = 5 
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Coupling constants: 

alpha.QED = 1 /( 137.036 ) 

alpha_s(Mz) = 0.1180 [ Mz = 91.1876 GeV ] 

alpha_s(iiiu) = 0.1709360043 

Quark masses : 

M_c = 1.65 GeV 

M_b = 4.75 GeV 

M_t = 175.00 GeV 

General switches (F=False, T=True) : 

only massless terms : F 

power suppressed terms : T 

QED corrections : T 

singlet contributions : T 

alphas'S m"2 included : T 

alphas~3 m"4 included : T 

alphas''4 m~2 included : T 
— end of parameters — 
R_had = 3 . 9049832 
R_u = 1.40765446 
R_d = 0.351913615 
R_s = 0.351913615 
R_c = 1.41577367 
R_b = 0.375554797 
R_t = . 
R_sing = -4.42977631E-05 
R_QED = 0.00221733842 

If Iverbose = .false., only the lines after " — end of parameters — " are printed. 
The parameter list displays the most important initializations, like the value of as (Mz) 
and the renormalization scale fj,. It also contains values for parameters that were derived 
from these settings, like the number of active flavors Uf, and the strong coupling constant 

10 Conclusions 

In this paper we have discussed the perturbative corrections to the inclusive cross section 
a(e'^e~ — > hadrons) and collected the analytical and semi-analytical formulas. This in- 
cludes the full mass dependence up to order a^, the expansion up to quartic mass terms 
at order and the quadratic mass corrections at order . Furthermore, the running and 
decoupling formalism necessary for a consistent evaluation of the strong coupling constant 
and the MS quark masses has been presented. 

The main subject of this paper, however, is a description of the f ortran program rhad 
which allows for the evaluation of R{s), including all currently available radiative correc- 
tions. The program is straightforward to use, in particular if the default parameter set 
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is adopted. Variations of the physical input values should be done such that they still 
resemble the physical case. The modularity of rhad allows for a simple extension once 
new corrections to the theoretical predictions for R{s) become available. In conclusion, 
rhad can easily be used to compute, for example, the perturbative parts of the hadronic 
contributions to the running of the electromagnetic coupling, and the anomalous magnetic 
moment of the muon. 
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Appendix 
A Installation 

The distribution of rhad contains the following files: 

Examples example. f makefile r012.f rhad.f vegas-rhad.f 

common. f funcs.f parameters. f r34.f runal.f 

Examples is a directory that contains programs to reproduce Tables 5, 6, and 7 of this 
paper. One example program is kept in the main directory, example. f. Its listing in 
shown in Sect. 9. It can be compiled by simply calling GNU make: 

> gmake prog=example 

The executable is named xexample. 

If GNU make is not available, the fortran files can be compiled individually using 
f77 -c -o <file>.o <file>.f 
and then linked using 
f77 -o xexample *.o 

B Basic functions of rhad 

Unless indicated otherwise, we use the following implicit type specifications in rhad: 
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IMPLICIT REAL*8(A-H,M-Z) 

IMPLICIT INTEGERd.J) 

IMPLICIT CHARACTER*60(K) 

IMPLICIT LOGICAL (L) 

The functions that encode the different loop orders to R{s) have already been discussed 
in Sects. 3-7, and in particular in Tables 1-4. 

The proper sums of these contributions are taken by the following functions corresponding 
to Rq{s), i?sing(s) and SR^^{s) (photon exchange only): 

rdqrk(s), ruqrk(s), rsqrk(s), rcqrk(s), rbqrk(s), rtqrk(s): 

Rd{s), Ru{s), Rs{s), Rcis), Rbis), Rtis) 

rsinglet(s): i?sing(s) 
rqed(s): Eq'^^§''°(«) 

The full i?-ratio is obtained by calling rhad(s), which adds the contributions from the 
individual quarks, depending on the c.m.s. energy. 

The general structure of a f ortran program from which the above functions are called is 
as follows: 

program <name of prograiii> 
<declaration of variables> 
<definition of s, e.g.:> 
include ' common. f 
scms = 11.5d0**2 
call parameters (scms) 
call init(scms) 

<call of function to compute R(s)> 
end 

It is important to call the subroutines parameters and init (in this order), which define 
the parameters and evaluate ctaip) and, if required, evolve and convert the masses. An 
explicit example program has been given in Sect. 9. 

C The subroutine parameters 
C.l Description of the parameters 

The file parameters. f contains the subroutine parameters and collects the most impor- 
tant variables that can be adjusted by the user, parameters. f should be understood as 
the input file for rhad. In Appendix C.2, our default settings are listed. Let us stress that, 
although many of the implemented formulas arc valid over a large parameter range, one 
should vary the input parameters only within a sensible region about their default values. 
Unreasonable settings, like inverted mass hierarchies or the like, may lead to inconsistent 
results. 
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The following list describes the variables defined in parameters. f: 

Iverbose: Print values of parameters (see Sect. 8) 
iunit: Output unit for parameter list. 

iord: Order of the calculation. Values: iord=0 (Born) to iord=4. It also governs the 
order at which the running coupling and the MS masses are evaluated. 

alphasmz: ai'\Mz)- Input value for the strong coupling constant (see Sect. 2). 
Imsbar: If .true., evaluate R{s) by using MS masses (see Sect. 2). 

masse, massb, masst: Initial values for mc, m;,, m^. If Imsbar = . true .: scale invariant 
mass ■fhq{rhq); otherwise: on-shell mass Mq (see Sect. 2). 

mu: Renormalization scale at which ag, the MS quark masses and R{s) are evaluated 
(see Sect. 2). 

muc,mub,mut: /ic, IJ'hi l^'t- Matching scales for the decoupling of the heavy quarks (see 
Sect. 2). 

Iqed: If .true., include the lowest order QED and, if iord.ge.2, the mixed QED/QCD 
corrections (see Eqs. (19) and (23)). 

Imassless: If .true., neglect all quark masses. 

Ipsup: If .true., include power-suppressed terms (see Eqs. (51), (54)). 

Ia3m2: If .true., include afmQ terms (see Sect.E.2). 

Ia3m4: If .true., include afmg terms (see Sect.E.2). 

laSsing: If .true., include singlet terms (see Sect. 6.2). 

Ia4m2: If .true., include a^rnq terms (see Sect. 7). 

thrc,thrb,thrt: s*^"", sf'"" (sec Sect. 1 and 8). 

Thresholds for open (perturbative) c-, 6-, and t-quark pair production. 

thrclow,thrblow,thrtlow: 4°"' 4°"^ (^ec Sect. 1 and 8). 

Perturbation theory is not applicable for Sq™ < s < Q G {c, 6, t}. Note that 

these three parameters have no influence on the results; their only task is to remind 
the user of the problematic energy regions in the perturbative approach. 

sqmin: minimal allowed value for the c.m.s. energy (see Sect. 1). 

Note that not all parameters are independent: for example, if Imassless = .true., 
Ia3m2, Ia3m4, and Ia4m2 will automatically be set to .false.. 
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C.2 Default settings 



The default settings in parameters. f are as follows: 



subroutine parameters (scms) 

c . . 

c. User-defined parameters, 
c . . 

implicit real*8(a-h,iii-z) 
implicit integer (i,j) 
implicit character*60(k) 
implicit logical (1) 



include ' common. f 



c . . verbose mode : 

Iverbose = .true. 



c. output unit for parameter list (6 = STDOUT) 
iimit = 6 



c. order or calculation: 
lord = 4 



c. strong coupling constant at scale mz (5 active flavors): 
alphasmz = O.llSdO 

c. use MS-bar or pole quark mass? (.true. == MS-bar mass) 
Imsbar = .false. 



masses 
masse 
massb 
mas St 



= 1.65d0 

= 4.75d0 
= 175. dO 



charm 

bottom 

top 



renormalization scale : 
mu = dsqrt(scms) 



decoupling scales: 

muc = 2.d0*massc ! charm 

mub = massb ! bottom 

mut = mas St ! top 



c. threshold for open quark production 
thrc = 4.8dO 
thrb = 11.2d0 
thrt = 2*masst+10.d0 
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c . . lower bound of quark threshold region 
thrclow = 3.73d0 
thrblow = 10.52d0 
thrtlow = 2*masst-10.d0 



minimum allowed cms energy: 




sqmin = 


l.SdO 




some switches 




Imassless = 


.false. 


use only massless approximation 


Iqed = 


.true . 


QED corrections (.true. = ON) 


Ipsup = 


.true . 


include power suppressed terms 


Ia3m2 


. true . 


\alpha_s~3 * m~2 terms 


Ia3m4 


.true . 


\alpha_s~3 * m~4 terms 


laSsing = 


.true . 


singlet contribution at order \alpha 


Ia4m2 


.true . 


\alpha_s~4 * m"2 terms 



return 
end 

D Implementation of running and decoupling 

Determination of OLa- rhad contains several routines wiiicli cover tlie consistent run- 
ning and decoupling of the strong coupling ai^^\^jL) (in this context see also Refs. [31, 18]). 

runalpha (apiO , muO ,mu , inf , inloop , verb , apiout ) : 

apiO: a^J^^\fio)/Tr 
muO: /xo 
mu: fi 
inf: Uf 

inloop: number of loops used for QCD f3 function 
verb: 0=quiet, l=verbose 
apiout: a^P^\^)/iT 

Evaluates a'J^^\ii) from ai^^\iiQ) by solving numerically the renormalization group 
equation 

decalphaCals .massth ,muth , inf , inloop , idir , alsout) : 
inf: Hi 

massth: M^, heavy quark mass to be decoupled 
muth: 

inloop: kord 
idir: 5 G {-1,1} 
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als: as ' 
alsout: as 



Evaluates as (/U**^) from aT {^J^^)■ The order of the matching relations is de- 
termined by fcord (0 < Aiord ^ 3). The matching coefficients require the pole quark 
mass as input. 

r unde calpha(alsO,muO,mul,alsout): 

muO: /Lto 
mul: /Lti 

alsO: ai"'^(/io)) where is defined in the subroutine parameters 

alsout: a^^^\ii{), where n/ is determined with the help of the variables thrc, thrb 
and thrt 

Determines a^^^\ixi) from Q;i"''^(/Lto), including decoupling and matching, by call- 
ing decalpha and runalpha. The value for rij =infini is set in the subroutines 
parameters, n/ =inf f in is determined according to the value of mul and becomes 
part of the common-block. 

In addition there is the subroutine decalphams which works analoguously to decalpha, 
except that the matching coefficients are parameterized in terms of MS quark masses. 

Determination of fhQ{y/s). 

runmass (massO , apiO , apif , inf , inloop .massout) : 

massO: rfiQ^fio) 
apiO: as(/Uo)/vr 
apif: as{i^f)/n 
inf: Uf 

inloop: number of loops 
massout: fhQ{jXf) 

Evaluates ffiQ^\ix) from •fnQ^\nQ) using Eq. (12). 
de cmas sms (mq , al s , mas sth , muth , inf , inloop , idir , mqout ) : 

mq: mQ"\fi^^) 
als: at\l^''') 

massth: rhh, MS value of heavy quark mass to be decoupled 
muth: 12^^^ decoupling scale 
inf: rii 

idir: S G {-1,1} 
inloop: kord 
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mqout: mJ^^'+'^V*'') 

Evaluates fhQ'''^^\ijL^^) from fnQ^\ij}^). The order of the matching relations is de- 
termined by kord (0 < kord < 3). The matching coefficients require the MS quark 
mass as input. 

rundecmass (mqO , inf ,muO ,mul , mqout) : 

mqO: rhQ'\iJ,o) 
inf: Hi 
muO: jUo 
mul: III 

_ (rif) / \ 

mqout: mg (/Hij 

Determines the MS mass fhQ^\fj,i) from ifiQ^\nQ), including decoupling and match- 
ing, by calling decmassms and runmass. n/ is determined according to the value of 
mul. 

Note that in decmassms and rundecmass only those cases are implemented which are 
needed for rhad. 

Pole mass versus MS mass. In case the MS quark masses are used as input the 
numerical values for the pole quark masses, which are needed for certain (double-bubble) 
contributions as described in Sect. 2, are computed with the help of the routine mms2mos. 

mms2mos (mms , inf , inloop ,mos) : 
mms: fnQ{fhQ) 

inf: number of active flavors 

inloop: number of loops used for the conversion 

mos: Mq 

Determines the on-shell mass Mq from the scale-invariant mass ffiQ^ffiQ), where the 
inloop-loop relation is used. In practice we use inloop=iord. For the conversion 
also asifhq) is needed which we determine with the help of the subroutines runalpha 
and decalphams. 



E Analytical results 

All results presented in this section are parametrized in terms of pole quark masses denoted 
by Mq, Ml or M2. The conversion to the MS formulas can be performed with the help of 
Eq. (67). 
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E.l Analytical results for Rq\s) 

In this appendix we list the results for the two-loop contributions to R{s) as they are 
implemented in rhad. They are either known analytically, or in terms of semi-analytical 
Pade approximants. 



The Abelian two-loop contribution reads 



(36) 



where r^^" is known in the form of Pade approximants [17]. We quote the result based 
on a [5/5] -Pade approximant: 



V — 



47r 



A8z 



(687 + 186^) Vl - + 216 In 



i-yF^T^ 



+ 



(1 + 
1- u 



1 + v/T^nyi 

2659.1447467995658619093787 + 994.90626341783737509012835 d) 

- 2742.9905230464882915790412 - 779.40067776480109262194206 
+ 562.98408428934139747541229 u)^ + 132.845153672444424204018458 



121.9768095064564342458129488 + 11.5590492918470301261688922 w 

- 119.511603603837393164213484 1^'^ - 1.04113850838776740702415214 lD^ 

-^ 18.4653468020861357442841578^"^ w^l }, (37) 



and 



ry^ can be found in Eq. (18). 
The non-Abelian contribution can be written in the form [17] 



z = 



1-vi 



'1- 



~ 2,^/1 1 
uj = l + 2i 



2 ■ 



5=4 



where ^ is the QCD gauge parameter, Vy^ is the result from the gluonic double-bubble 
diagram given below and r^'^" is again given in terms of a Pade approximant: 



(38) 



(39) 
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A,pa 

V = 



CaCf — Im 
47r 



33 - 27 z - 6^2 ^ 

ZTT ; , + [1+ UJ 



8zVi- lA 

53.56935816521245 + 23.42800818317322 - 28.58891755403711 



10.11103466811403 w''^ + 0.4266396586616486 lD^ 



19.19224589510177 - 1.767625499084376 w - 13.2990319119615 

i1 



+ 1.00664474828887 cD^ + 



(40) 



Following Ref. [32], the exact results of the imaginary part of the fermionic and gluonic 
double-bubble diagram (see Fig. 3 (d) and (e)), ry and Vy'^, respectively, may be written 
as (x G {L; A,^}) 



with 



- — ( RZ^ It^— — Rn 



+ RZ.S('^ 



(41) 



Cl = 

CA,g = 

rIu = 



J. It-^ 



R. 



CfT, 
CaCy , 

1, 

-^ + ln4, 

5 3 

4 8^' 

12 4^ 16^ V 4 8^ ' 



(42) 



Vy^ is defined in (18), whereas S^^^ was originally calculated in Ref. [7] and reads: 
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c(2) ^ (3-V) 

6 

X |Li3(p) - 2Li3(l -p) - 3Li3(p2) - 4Li3(^) - 5Li3(l -p^) + 11^3 



+ Li2 (p) In 



4(1-V)- 



■1- V 



+-^m 



361n21np- 441n^p + 491npln (^-^^ 
36 In 2 + 21 Inp + 16 In vq - 22 ln(l - vq"^) 
+ ^|(15 - 6vQ^ - VQ*) (hhip) + Li2(p2)) + 3(7 - 22vq' + 7vQ*)U2{p) 



AvQ 
1- V 



-lnplni;Q 



+ 



(1 - VQ){bl - 45t;Q - 21vQ^ + 5t;Q3)C2 

(1 + vq) (-9 + 33i;q - 9^Q^ - 15i;q=^ + Avq-^) 
VQ 



+ 



(33 + 22vQ^ - IvQ^) In 2 - 10(3 - vq^){\ + vq^) Xt^vq 



(15 - 22i;q2 + 3^q4) 



1- V 

4wq2 



Inp 



+ 2i;Q(3- 



4(1-V) 

,2 I oo„,^3 



^ 237 - 96.^ + 62 V + 32V -59V _ ,6,^(3 _ ,^2^ ( 1+^ 

i;Q(75- 29V) 



2.;Q(39-17V)ln(^^^)- 



(43) 



with C2 = 7rV6 ~ 1.64493, Cs ~ 1.20206, and 



1 + VQ 



The double-bubble contribution Vy can be written in the form 

/ 1 Ml r^^) ^ 

tI{s,Mq) = TCf I py{M^,M^,s) +p^{Mi^,M^,s) + -In 



where 



P^{M'q,M^,s) 



3 + lOf 2 - „ 
^ + 



-3 + 40i;| + lQv% - lbv% 



Q 



+ 



+ 



-18 + 234w^ + 167u^ - 118u| 



124 
-3 - IQvi, + ^vi 



In p 



(44) 



-9 + SlOv^ - 118?;^ 



+ 



Q 



+ VQ {-27 + bvl)C2 



(45) 
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and 



P^{MlMls) = \ [ 
with 



(1-2MQ/Vi)" 
4M2/S 



r(i-Vy? dz ( 2Ml 
dy — 1 + 

JiMl/s Z y SZ 



4M| 
sz 



(46) 



8M4 



^{y, z) 



^ + ^{l-y + z)-{l-y + zf-2{l + z)y 
1-y + z 



X In ■ 



1-y+z- 



sy 



4M2 



l-y + 2:+^/l-^AV2(l,y,^) 



A(l,y,2) 



^A'/\l,y,z) 

V 



1 + y2 + - 2(y + z + yz) . 



16M; 



Q I 
^ r 



+ 4 1 + 



2M2 



(l-y + z)2 



4M^ 



A(l,y,^;) 



(47) 



Note that p^{M^,M^,s) vanishes for s < (4Mq)2. 

The numerical integration of Eq. (46) would be a bottle neck of rhad. But it turns out that, 
for s > (4Mq)^, (s, Mq) is extremely well approximated by the high energy expansion of 
Ref. [16], and it drastically shortens the running time of rhad. Thus we use this expression 
as the default in rhad for s > (4Mq)^, whereas for s < (4Mq)^ we can savely use Eq. (44) 
since the double integral does not contribute in this case. For completeness, we list the 
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expansion for ry up to Mq/s° [16]: 
4{s, Mq) = Cft| - y + Cs + \Ls^ + ^ 



^7 + ^Lsn 



+ 



+ 



+ 



+ 



+ 



- + I8C2 - IOC3 + 12L„. - 3L^, + U - ^^ms L 



4 352 

77 + ^ 



350 



76 



188 116 



ms I -^s/i 
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20233 533 1673 _ 25 2 
^88" + "6" ^ ~W ~ T 



983 203 



36 



54559 3592 



(48) 
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1157 



4328 



+ - 



6750 45 
61699 4676 



-Li 



675 



45 



ms j ^sfj. 



'Ml' 



+ 



9214697 , 346981 18937 . 
1105 C2 + -^Lms + ^g-^L 



6480 
84743 3064 



270 9 
2 /„\ r _ /,,2 



+ 



where = ln(MQ/s), Lg^ = ln(s//i ). The ellipse denotes uncalculated higher order 
terms in Mq/s. 

The double-bubble contribution with zero outer mass and inner mass M2 is given by [9] 



rf(s,0,M2) = CFT 



q"" {Ml s) + 6{s - 4M|) ( Q^'iMl s) + -\u 



1 . M?: 



(49) 
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with 



q\mIs) = ^(1-6x2) 



lLi3(l^)-lLi3(i±^: 
2^2' 2^2' 



l + w 1-w l + w l-w 



2 'l-w' [ 12 'l-w' 2 a + 1' 

-lln(i±^)ln(i^' 
2 ^ 2 2 ' 

1 / . x/^./l + 'fi'x ^. A — w^ ^ A + w. ^. ,1 — w. 
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''a + 1 1 — w 
19. 2\ f . r 1 + w^ ^. , 1 — w^ , ,,l + i(;. 



+ i[—+x + x') Li2(-- )-Li2(-— — )-lnxln( 

\72 J \ 1 — w l + w 1 — w 

^/ 73 74 2^,/l + «^^ 1/2123 2489 \ 

/(M|, s) = ^ (1 - 6x2) J^Li3(^=^) - Cs - 2C2 In^ + ^ In^ 

+ ^ (19 + 46x) ^/l+4^ (Li2(^2) - C2 + In^ A) 
5/53 ^ \ , 3355 119 



where 



X = , a = vT+4x, ■u; = Vl - 4a; , A= "^^^"^^ ^ 



s V4x 



In the hmit s » M|, the result for ry'(s,0,M2) reads 



11 . 1, s „ /Ml 
-In- 

4 (Lt 



r*(.,0,M2) = CFT(-- + C3 + 7lnT2+^^) ) ■ (^2) 



Furthermore we need ry'{s, Mi,0) up to the quadratic mass terms which enters the eval- 
uation of R^^ in the limit s > sf^. It is given by 



(53) 

In order to obtain the result for the double-bubble contribution in the limit M| ^ s one 
has to perform the decoupling of M2 from the coupling constant using Eq. (9) and the 
one-loop approximation of Eq. (65) . Taking into account the full Mi dependence the result 
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can be cast in the form [11] 



4'^'\s,Mi,M2) = CfT 



l2Mf 16Mf\, l + A/l 
■1 + ^ + ^ 1 In- 



4Mf 



l-\ 1 



AMI 



(22 79M2 SMf\ r 



4M2 



+ 



(54) 



fl + ^Vn4\/^ 
V s J s 

The expansion of CpTg^ in the hmit S> s agrees with r^'^ for Mi = and reads 

db,V HM^ n ^ 8 M|\ 

r^'(.,0,M,) = ^^FT^I^^ + ^ln— J . (55) 

The massless hmits of the above expressions have been given in Eq. (22) . 



E.2 Analytical results for Rq\s) 

This section hsts analytical expressions for the three-loop functions defined in Sect. 6. 
As already mentioned, only the high-energy expansion is known, including quartic mass 
terms. We write these functions as follows: 



,(3) _ .(3) 



'0 



+ 



J3) 
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,(3) 

y,4 ' 
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„(3) 



(56) 



The results for the massless terms have been obtained in Ref. [21, 22], the quadratic mass 
terms in Ref. [23] , and the quartic terms in Ref. [24] . Adopting the on-shell scheme for the 
mass Mq, we get, for the contributions where the massive quark couples to the external 
current: 
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(59) 



where n/ is the number of active flavors, 04 = Li4(l/2) « 0.517479, C2 = vr^/6 « 1.64492, 
Ca 1.20206, C4 = 7r'^/90 ^ 1.08232, Cs » 1.03693, L^^ = ln(M|/s) and = ln(s///2). 
The mass corrections in the case where the massive quark only appears as insertion in 
gluon lines reads: 
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(61) 



E.3 Analytical results for Rq\s) 



With the help of the rcnormalization group equation the logarithmic contributions of the 
massless order result can be reconstructed exactly. Thus, for fJ,^ ^ s Eqs. (33) and (34) 
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take the form 
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where the first number in r, 



(4),0L (4),1L -(4),0L 
' '0 



and ry2^^ corresponds to the estimates 



obtained with the help of the results of Refs. [26, 27]. 



E.4 Renormalization group coefficients 

This section collects the formulas used for the evaluation of the running coupling constant 
and the MS quark masses. It also contains the conversion formulas for transforming MS 
quark masses to their on-shell values. 

The coefficients of the /3-function defined in Eq. (8) read [33] : 
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The coefficients of the unction of Eq. (10) are [34]: 
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The decoupling coefficient for of Eq. (9) is [35]: 
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and for the quark masses, Eq. (13) [35]: 
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(66) 



where S4 ~ —1.76280, n/ = — 1 is the number of light (massless) quarks and M^ is the 
pole mass of the heavy quark. The versions of Eqs. (65) and (66) where the MS mass is 
used as parameter can be found in Refs. [36, 35, 31, 18]. 

The conversion from MS to on-shell quark masses reads [3, 4, 5]: 
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where 04 = Li4(l/2) 0.517479, l^m = ln(/x2/mQ(//)), and ri; = ra/ — 1. Other variants of 
Eq. (67) are listed in Ref. [31]. 
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